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I. INTRODUCTION 


Electromagnetic interactions have long been used to study both hadronic structure and 
strong interaction dynamics. Examples include deeply inelastic lepton-nucleon scattering, 
hadropro duct ion of lepton pairs, the production of photons with large transverse momenta, 
and various photoproduction processes involving scattering of real or very low mass virtual 
photons from hadrons. In particular, heavy quark production in deeply inelastic electron- 
proton scattering (DIS) is calculable in QCD and provides information on the gluonic content 
of the proton which is complementary to that obtained in direct photon production or 
structure function scaling violations. In addition, it forces one to address such issues as 
when the mass of the heavy quark should be neglected, and how this is done consistently. 

Activity in the area of neutral current DIS charm quark production has increased re¬ 
cently with new data becoming available from ZEUS [Q and HI p| at HERA. In particular, 
substantial samples of D*^(2010) hadrons have been obtained. On the theoretical side, 
much attention has been given to the study of heavy quark contributions to the proton 
structure function. They have been calculated to next-to-leading order in fully inclusive 
0, single inclusive Q], and fully differential (exclusive) 0 forms as QCD corrections to 
the photon-gluon fusion process (i.e. in -|- 0{a^) using three flavor parton densities). 

Experimentally and phenomenologically [|^Q, it is this process that is seen to dom¬ 

inate near the threshold region. However, well above threshold the heavy quark may be 
considered massless and included in the parton distribution function of the proton. Various 
schemes to match these two regions have been proposed MI. Despite all of the attention 
structure functions have received, much less has been done to explore the actual heavy quark 
differential cross sections in DIS which are, in fact, much easier to measure experimentally. 

The purpose of this work is to present next-to-leading order cross sections for charm quark 
production in the x and region covered by the HERA collider. Additionally, predictions 
for heavy hadrons, namely D*^(2010), are given. The calculation is implemented in a Monte 
Carlo style program which allows the simultaneous histogramming of many distributions 
incorporating experimental cuts. It represents an elaboration of the brief results already 
presented by one of us [|^, and an application of the fully differential heavy quark structure 
functions calculated in [^. No experimental data is shown here because the results have 
already been added to several of the plots in (see also 0)- An extensive comparison will 
be given elsewhere. Herein, the calculation itself is discussed and the variation with respect 
to the theoretical parameters is studied. Heavy quark correlations have also been calculated 
for hadropro duct ion |^, photoproduction |^, and photon-photon collisions [|^,^ allowing 
for the possiblility of an extensive comparison with experimental data. 

The calculation was performed using the subtraction method which is based on the 
replacement of divergent (soft or collinear) terms in the squared matrix elements by general¬ 
ized distributions. The method was first used in the context of electron-positron annihilation 


T7| and its essence is described and compared to the phase-space slicing method in the 
introduction of a paper by Kunszt and Soper M . 


The remainder of the paper is as follows. A review of the subtraction method and other 
aspects of the calculation, including how the hadronization is modeled, are given in Sec. H. 
Numerical results and a discussion of related physics issues are presented in Sec. HI. The 
conclusions are given in Sec. IV. 
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II. CALCULATION 


In this section the calculation of the charm quark cross section in deeply inelastic scatter¬ 
ing is described. First, the cross section is written in terms of the charm quark contribution 
to the proton structure functions. Then the next-to-leading order QCD corrections to the 
structure functions are reviewed. Finally we describe the connection with the production of 
heavy hadrons containing a charm quark. 


A. Cross section in terms of strncture fnnctions 

The reaction under consideration is charm quark production via neutral-current electron- 
proton scattering. 


e (/) -1- P{p) —> e (I') + Qipi) + X . 


( 2 . 1 ) 


When the momentum transfer squared = —q^ > 0 (g = / — /') is not too large -C M|, 
the contribution from Z boson exchange is kinematically suppressed and the process is 
dominated by virtual-photon exchange. The cross section may then be written in terms 
of structure functions Q^,m) and F£(x,Q^,m) which depend explicitly on the charm 
quark mass rric | 

(fa 


as follows: 
27rQ:^ 


dydQ"^ yQ^ 


{[l + (1 - v)'^] Q'\ “c) - V^FHx, Q^, m^)} 


( 2 . 2 ) 


where x = 1‘ly ■ q and y = p ■ q/p ■ I are the usual Bjorken scaling variables and a is the 
electromagnetic coupling. The scaling variables are related to the square of the c.m. energy 
of the electron-proton system S = (/ + via xyS = Q^. The total cross section 13 is 
given by 


JAraljS Jraly^l{\-y) XdydQ^ 


(2.3) 


where rUe is the electron mass. In deriving Eq. (|2.2|) , one integrates over the azimuthal angle 
between the plane containing the incoming and outgoing electrons and the plane containing 
the incoming proton and the outgoing charm quark. 

As mentioned in the introduction, this process is described near threshold in the frame¬ 
work of perturbative QCD by flavour creation through the virtual-photon-gluon fusion pro¬ 
cess 


Y{q) + gih) ^ Q(,Pi)+Q{ p 2 )- (2.4) 

The structure functions follow from the longitudinal a^ and transverse aT cross sections 
for this reaction via = (Q^/dvr^a) (<7^ -|- ax) and = (Q^/dvr^a) Thus, QCD 
corrections to the reaction (0 correspond to QCD corrections to (|2.4|) to which we now 
turn. The supersript c will henceforth be dropped to simplify notation. 
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B. QCD corrections to the heavy quark structure functions 


Within the context of perturbative QCD, structure functions are expressed as a product 
of the running coupling, the parton densities, and the hard scattering cross sections. The 
result is a physical quantity, but the individual terms can be defined in a convenient scheme 
which moves terms from one factor to another. All schemes should give the same result for 
the product, up to terms of higher order. 

A next-to-leading order calculation of the heavy quark contributions to the proton struc¬ 


ture functions requires the one-loop virtual corrections to (|2.4|) . For this set of diagrams, 
the renormalization was carried out so that divergences coming from the light quarks were 
subtracted in the standard MS scheme, while the divergences coming from heavy quark loops 
were subtracted at zero external momenta. This is the scheme originally proposed in 


in which the mass m only appears in the hard scattering cross sections. As a consequence, 
below the subtraction scale, only the number of (massless) light flavours appears in the 
running coupling and in the splitting functions used in the parton evolution equations. At 
the subtraction scale fi = m there are matching conditions for both the running coupling 
and the light flavour densities Therefore, to order there is no charm density at the 


scale = m B. [0 


Consequently, one should use only parton distribution sets that were derived from data 
using the same renormalization scheme. Examples of such densities are GRV94 and 
CTEQ4F3 At larger scales there is a charm density proportional to as ln(/i^/m^), which 
grows at the expense of a reduction in the gluon density. One of the interesting problems in 
the analysis of charm quark contributions to deeply inelastic scattering is to understand the 
transition region from the photon-gluon predictions based on three light flavours to a charm 
density picture with four light flavours. The matching conditions become more complicated 
as one goes to higher order. The corresponding two-loop matching conditions for the flavour 
densities have been calculated in [^] wherein they were found to have finite terms at fi = m. 

In addition to the virtual corrections described above, there is also a contribution from 
the gluon-bremsstrahlung process 


7'(i) + ^ Q(pi) + Q(P 2 ) + s(h) 


(2.6) 


and new production mechanisms not present at leading order, which are given by 

j*(g) + g(ki) Q(pi) + Q(p 2 ) + g(k 2 ) 

j*(g) + g(ki) Q(pi) + Q(p 2 ) + ^(^ 2 ) (2.6) 

where (g)g is a massless (anti-)quark. The contribution to the structure functions resulting 
from these processes have been calculated in a fully differential P| form and are suitable 
for use in constructing a Monte Carlo style program because one has full access to the final 
state partonic four vectors. 

Briefly, the computation in was carried out using the subtraction method which 
is based on the replacement of divergent (soft or collinear) terms in the squared matrix 
elements by generalized plus distributions. This allows one to isolate the soft and collinear 
poles within the framework of dimensional regularization without calculating all the phase 
space integrals in a space-time dimension n 7 ^ 4 as is required in a traditional single particle 
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inclusive computation. In this method the expressions for the squared matrix elements in the 
collinear limit appear in a factorized form, where poles in n — 4 multiply splitting functions 
and lower order squared matrix elements. The cancellation of collinear singularities is then 
performed using mass factorization. The expressions for the squared matrix elements in 
the soft limit also appear in a factorized form where poles in n — 4 multiply lower order 
squared matrix elements. The cancellation of soft singularities takes place upon adding 
the contributions from the renormalized virtual diagrams. Since the final result is in four¬ 
dimensional space time, one can compute all relevant phase space integrations using standard 
Monte Carlo integration techniques. 

The resultant (differential) structure functions may be written in the form 


Fk{x,Q‘^,m) = 








ctl + 47ra,(/i^) ( In 


( 1 ) 










i=q,q 


ei ( c2 + eg In ^ 1 + ef dl"} + en o)!} 


2 Ti) 






(2.7) 


k = 2, L. The lower boundary on the integration is ^rain = + Q‘^)/Q‘^. The parton 

momentum distributions in the proton are denoted by fi/p{^, The mass factorization 
scale /i/ has been set equal to the renormalization scale and is denoted by /i. All charges 


are in units of e. Finally, eg, c^, , (i = g,q,q), and d^^j, = q,q) are scale 

independent parton coefficient functions. They are functions of and m. In Eq. 

( p.7|) the coefficient functions are distinguished by their origin. The c-coefficent functions 
originate from processes involving the virtual photon-heavy quark coupling, while the d- 
coefficient functions arise from processes involving the virtual photon-light quark coupling 
and the o-coefficent functions are from the interference between these processes. 

In addition to the calculation of [^, the functions and were calculated in 

inclusive form in PI, [01 as two-dimensional integrals and computed numerically. In E3 


(b ^(b 


they were numerically tabulated in grids, with a fast interpolation routine, so that the 
computation of Eq. (|2.71) could be included in a global fit, if desired. Finally, in p6 


exact 


analytic formulae were given for the d^^\ together with analytic formulae for all the coefficient 


functions c 


( 0 ) 
k.i ’ 


db 


db 

''k,i 


i{'i' = 9 i Qj 9)1 and ,{i = q, q) in the limit Q 


^ >> w?. 


The latter 


results are necessary to consider a variable flavour scheme in which the coefficient functions 
are incorporated into redefined light parton densities including a charm density [|n . 


Naturally, properties of the structure functions give insight into the behavior of the 
cross section. Therefore the more salient features of the next-to-leading order structure 
functions will now be summarized. The interested reader can find additional details in the 


P7H21 


For 
10-^ is 


original papers and, more so, in the recent phenomenological analyses 
moderate ~ 10 GeV^ one finds that the charm quark contribution at small x 
approximately 25% of the total structure function (defined as light parton plus heavy quark 
contributions). In contrast, the contribution from bottom quarks is only a few percent 
due to charge and phase space suppression. Thus, the charm quark contribution must 
be retained, but the bottom quark contribution may be neglected in the following. The 
gluon initiated contributions (p2.4|) and (|2.5|) comprise most of the structure function. The 
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FIG. 1. F^^x, Q‘^,mc = 1-5 GeV) as a function of x for GRV94 (solid lines) and CTEQ4F3 (dash 
lines) parton distribution sets for = 3GeV^ (bottom), 25GeV^ (middle), and 50GeV^ (top). 


quark initiated processes (|2.6|) give only a few percent contribution at small x for reasonable 
scale choices. Results for the charm quark contribution to the proton structure function 
F 2 ^(x, Q^,mc = 1.5 GeV) as a function of x are shown in FIG. 1 for /i = \JQ‘^ + 4m^ using 
the GRV94 |2^ and GTEQ4F3 iQ parton distribution sets. As mentioned above, these 


sets were chosen because they have rij = 3 in the evolution and are therefore the most 
consistent sets to use with the NLO calculation. The curves show a marked rise in the 
structure function at small x due primarily to the rapidly rising gluon distribution. The 
renormalization/factorization scale dependence is rather flat, especially at small x where 
the structure function is largest. This is demonstrated in FIG. 2 for various x and 
values. By far the largest uncertainty in the calculation of the structure functions is the 
value of the heavy quark mass. For charm production, for example, a ±10% variation of 
the mass about the central value of rric = 1.5 GeV gives a ±20% variation in the structure 
function for moderate Q^. 


C. Fragmentation into heavy hadrons 


Experimentally, it is heavy hadrons or their decay productions that are observed. Of 
interest for HERA is 11*^(2010) production. Herein, the Peterson et al. form of the frag¬ 
mentation function [^] is used to model the nonperturbative hadronization process. 

The cross section for D* production is obtained by convoluting the charm quark cross 
section (p^.3D with the fragmentation function 


D(z) = 


N 




( 2 , 8 ) 
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II (GeV) II (GeV) 

FIG. 2. ,mc = 1.5 GeV) as a function of the scale /x with rric < /U < 2-\/ Q‘^ + 4m^ for 

various x and values for the GRV94 (solid lines) and CTEQ4F3 (dash lines) parton distribution 
sets. 


where N is fixed such that D{z) is normalized to unity. The normalization of the cross 
section is then fixed by the charm quark fragmentation probability P{c —>■ D*) = 0.26 


The parameter e may be extracted from data and used as input. However, in light of 
recent work on the subject [^,Q, the specihc value that should be used in this particular 
application is not obvious. The choice of the best value is left as the subject for future study. 
Below e is treated as a free parameter, and the effect of its variation on the cross section 
will be examined and considered as part of the uncertainty due to hadronization. 

Other sources of uncertainty include such technical details as how exactly the four vector 
of the D* is formed. One may scale the entire four vector by z, but then the hadron mass 
is zrric- Another possibility is to scale the three vector by z and £x the energy component 
such that that the mass is m{D*) = 2.01 GeV. The latter is used here. 

Evolution of the fragmentation function, which one expects to become important when 
Pt 3> rric, is not included because the region of interest is pt ~ rric- Indeed, recent calculations 
of charm photoproduction at HERA [^] have shown that this effect becomes sizeable only 
for Pt > 20 GeV. 


III. RESULTS 

Using the results of the method described in the previous section, a computer program 
has been constructed to calculate charm quark and/or D* cross sections in deeply inelastic 
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FIG. 3. Comparison of leading order results (open circles) with AROMA ||^ (solid lines) for 
charm production in the kinematic region 10 < < 100 GeV^ and 0 < y < 0.7 at HERA. 

scattering Q The program uses Monte Carlo integration so it is possible to study a variety 
of distributions and implement experimental cuts, provided they are dehned in terms of 
partonic variables or the optional fragmentation into heavy hadrons is used. Results are 
presented in the HERA lab frame with positive rapidity in the proton direction. The proton 
and electron beam energies are taken to be 820 GeV and 27.6 GeV, respectively. There 
are several necessary cross checks that should be performed before discussing the general 
properties of the complete next-to-leading order calculation. 


A. Comparison with existing results 


A comparison with the leading order event generator AROMA provides a check 
on the kinematics and leading order matrix elements through the shape and overall nor¬ 
malization of the distributions. For this purpose, both codes were run with the CTEQ2L 
3^ proton parton distributions, the default set for AROMA. This set has = 190 MeV 
which was used, along with uj = 4, in the one-loop strong coupling ag- The renormal¬ 
ization and factorization scales were both set to '/§. In the AROMA calling program all 
fragmentation and showering was turned off, and only the contribution from virtual photon 
exchange was retained. Shown in FIG. 3 are the results for charm quark production, as¬ 
suming rric = 1.5 GeV, in the kinematic region 10 < < 100 GeV^ and 0 < y < 0.7. The 

distributions shown are for the transverse momentum pt and pseudo-rapidity r] of the charm 


^Interested readers should contact harris@hep.fsu.edu for a copy of the computer code. 














FIG. 4. Leading (dash) and full next-to-leading (solid) order differential cross sections for charm 
quark production at ^/S = 301 GeV at HERA using the GRV94 (LO and HO, respectively) parton 
distribution set at /i = + 4m^ with nic = 1.5 GeV. 


quark, both in the HERA lab frame, along with the usual deep inelastic scattering variables 
X and Q^. These will be the canonical set of observables used in the rest of the paper. The 
shapes of the curves are virtually identical over several orders of magnitude and the area 
under the curves is the same to better than 1%. 

Another check is to reproduce the total heavy quark cross section as previously calculated 
in next-to-leading order |^. Both computer programs were run with the CTEQ3M 


proton parton distributions. This set has A^^^ = 239 MeV which was used, along with 
nj = 4, in the two-loop strong coupling The renormalization and factorization scales 
were both set to + 4m^. The results for charm quark production, again assuming 


rric = 1.5 GeV, in the kinematic region 10 < < 100 GeV^ and 0.01 < y < 0.7 are identical 

for both programs to better than three signihcant hgures. The numerical values for the cross 
section per channel are 7.48 nb for 0(as) photon-gluon, 2.68 nb for 0(a‘^) photon-gluon, and 
—0.41 nb for the sum of photon-quark and photon-antiquark. Having made these checks, 
the general properties of the full next-to-leading order cross section may now be considered. 


B. Properties of the charm quark cross section 


In this section the dependence of the charm quark cross sections will be studied as a 
function of parton distribution set, charm quark mass, and renormalization/factorization 
scale. All results are for the kinematic range 3 < < 50 GeV^ and 0.1 < ?/ < 0.7. 

The GTEQ4F3 and GRV94 HO proton-parton distribution sets were used in 
the remainder of the paper, as noted. For leading order results, the GRV94 LO p3 set 


9 































FIG. 5. Next-to-leading order differential cross sections for charm quark production at 
\/5 = 301 GeV at HERA using the GRV94 (dash) and CTEQ4F3 (solid) parton distribution 
sets at /r = y/Q‘^ + 4m^ with nic = 1.5 GeV. Also shown for comparison in part (c) is da/dlog{^) 
vs. log(^) (right pair of curves). 

was used. The (one-) two-loop version of the strong coupling was used with matching 
across quark thresholds for the (LO) NLO results. The value of Aqcd was taken from the 
proton-parton distribution set. The renormalization and factorization scales have been set 
equal to /r. 

The leading (dash) and next-to-leading (solid) order differential cross sections for charm 
quark production using the GRV94 (LO and HO, respectively) parton distribution set at 
/i = \JQ'^ + 4m^ with rUc = 1.5 GeV are shown in FIG. 4. The shape of the NLO transverse 
momentum distribution is similar to the LO one, but somewhat flatter. The pseudo-rapidity 
distribution shows the radiative corrections are concentrated in the negative rapidity direc¬ 
tion, tending to pull the maximum back towards the central region. The Bjorken x distri¬ 
bution receives corrections near its maximum with near zero correction at the tails. The 

distribution receives a fairly uniform shift in normalization. Save at high pt, the NLO 
predictions lie below the LO ones. This a property of the GRV parton distribution set. 
For the GTEQ4F3 set the opposite behavior is observed. This a reflection of the difference 
between the leading order gluon distribution functions and the corresponding Aqcd of the 
two sets. 

One may ask how sensitive are the full next-to-leading order results to modern parton 
distribution sets. The answer is immediate from FIG. 5 which shows the next-to-leading 
order differential cross sections for charm quark production using the GRV94 (dash) and 
GTEQ4F3 (solid) parton distribution sets at /i = with irtc = 1.5 GeV. From Eq. 

( p.7|) the parton distributions are probed at ^ which is typically one order of magnitude larger 
the X. This is illustrated in FIG. 5c where a plot of da/d\og{^) vs. log(^) (right set of curves) 
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FIG. 6. Next-to-leading order differential cross sections for charm quark production at 
y/S = 301 GeV at HERA using the GTEQ4F3 parton distribution set at // = ^/Q'^ + 4m^ with 
rric = 1.35 GeV (solid) and rric = 1.65 GeV (dash). 

is superimposed on the plot of da/d\og{x) vs. log(a;) (left set of curves). The difference 
between the curves produced using the two parton distribution sets is approximately 10% 
at .^ = Thus, the predictions are less dependent on the parton density sets in NLO. 

The largest uncertainty in the structure function calculation is that of the charm quark 
mass. The same is true for the cross section as shown in FIG. 6 for the next-to-leading 
order differential cross sections for charm quark production using the CTEQ4F3 parton 
distribution set at p = \JQ‘^ + 4m^ with rric = 1.35 GeV (solid) and rric = 1.65 GeV (dash). 
Mass effects are smaller at the larger transverse mass because they are suppressed by powers 
of mclvt in the matrix elements. However, as mentioned above, if the range is extended much 
further, large logarithms of the form ln(p^/m^) appear in the cross section and should be 
resummed. 

Finally, the scale dependence is shown in FIG. 7. The next-to-leading order differential 
cross sections are for the GTEQ4F3 parton distribution set at p = 2mc (solid) and p = 
2^Q‘^ + 4m^ (dash) with rric = 1-5 GeV. The curves show very little scale dependence. This 
can be anticipated from the results shown in FIG. 2 and the distribntion in Bjorken x shown 
in FIG. 7c. The latter shows the cross section is dominated by x 10-3-® = 3.2 X 10-^ while 
the former shows that, independent of the structure function is very flat in this particular 
X region. Therefore, the cross section tends to be fairly insensitive to the choice of scale. 
Other kinematic regions show increased scale dependence. This serves as a reminder that 
care must be taken in interpreting the results of varying the renormalization/factorization 
scale to estimate the size of the theoretical error. 
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FIG. 7. Next-to-leading order differential cross sections for charm quark production at 
\/5 = 301 GeV at HERA using the CTEQ4F3 parton distribution set at /r = 2mc (solid) and 
^ = 2yjQ‘^ + 4m^ (dash) with rric = 1.5 GeV. 


C. Predictions for D* production 

In this last section the fragmentation is tnrned on and predictions for D*^ production at 
HERA are given. The power of the subtraction method becomes manifest because experi¬ 
mental cuts can easily be implemented. Cuts similar to those preferred by ZEUS [Q and HI 
0 are used. Namely, 2 < < 100 GeV^, 0.05 < y < 0.7, pf* > 1.5 GeV, and \t]^*\ < 1.5. 

No distinction is made between D*^ and D *~. The results shown in FIG. 8 use the GRV94 
parton distribution set at p = with nic = 1.5 GeV and e = 0.03 (dot), e = 0.06 

(solid), e = 0.09 (dash). The variation in the area under the curves is roughly half that 
from the mass uncertainty. The shape changes are very mild. Using the GTEQ4F3 parton 
distribution set instead would give slightly larger results. 


IV. CONCLUSION 

The calculation of next-to-leading order corrections has allowed more reliable predictions 
of heavy quark differential distributions in deeply inelastic scattering. In addition, with the 
calculational formalism used here, it was possible to add Peterson fragmentation thus giving 
predictions for 17*^(2010) production at HERA. 

At leading order the results were cross checked against AROMA and found to give 
complete agreement. When the program is run in fully inclusive mode it reproduces existing 
results for the charm quark cross section at next-to-leading order. 

The radiative corrections to the lowest order photon-gluon fusion process are important 
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log(x) (GeV=) 

FIG. 8. Next-to-leading order differential cross sections for D* production at ^/S = 301 GeV 
at HERA in the kinematic region 0.05 < y < 0.7, 2 < < 100 GeV^, 1.5 < Pf* < 10 GeV, and 
\r]^ I < 1.5 using the GRV94 parton distribution set at jj, = \/Q'^ + 4m^ with rric = 1.5 GeV and 
e = 0.03 (dot), e = 0.06 (solid), e = 0.09 (dash). 

as they change both the shape and normalization of the transverse momentum, pseudo¬ 
rapidity, and X distributions. The distributions only receive a shift in normalization. In 
the kinematic region studied, the cross section is very stable with respect to variations in the 
renormalization/factorization scale because the cross section is dominated by an x region 
where the scale dependence of the underlying structure functions is nearly flat. The scale 
dependence of the hard scattering cross sections is well compensated by that of the parton 
densities and Ug- 

The cross section is dominated by the rapidly growing gluon distribution at small x, but 
distinguishing between modern parton distribution sets using this process will be difficult as 
demonstrated by the fact that they give nearly identical results for a variety of observables. 
This is compounded by relatively large uncertainties from the quark mass and hadronization 
effects. At present a comparison with experimental data can offer a confirmation of the gluon 
distribution at small x. Examining a variety of x and bins will shed light on the transition 
region between massive and massless charm quark descriptions. 
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